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We have investigated dynamical heterogeneity in a highly supercooled liquid using molecular- 
dynamics simulations in three dimensions. Dynamical heterogeneity can be characterized by three 
quantities: correlation length £4, intensity X4, and lifetime Thetero- We evaluated all three quantities 
consistently from a single order parameter. In a previous study (H. Mizuno and R. Yamamoto, 
Phys. Rev. E 82, 030501(R) (2010)), we examined the lifetime Thetero(t) in two time intervals 
t — r a and T ngp , where r a is the a-relaxation time and T ngp is the time at which the non-Gaussian 
parameter of the Van Hove self-correlation function is maximized. In the present study, in addition 
to the lifetime Thetero (t), we evaluated the correlation length and the intensity Xi(t) from the 
same order parameter used for the lifetime Thetero (t). We found that as the temperature decreases, 
the lifetime Thetero (t) grows dramatically, whereas the correlation length £((£) and the intensity 
Xi{t) increase slowly compared to Thetero (t) or plateaus. Furthermore, we investigated the lifetime 
Thetero (t) in more detail. We examined the time-interval dependence of the lifetime Thetero (i) and 
found that as the time interval t increases, Thetero (i) monotonically becomes longer and plateaus at 
the relaxation time of the two-point density correlation function. At the large time intervals for 
which Thetero (t) plateaus, the heterogeneous dynamics migrate in space with a diffusion mechanism, 
such as the particle density. 

PACS numbers: 64.70.P-, 61.20.Lc, 61.43.Fs 



I. INTRODUCTION 

As liquids are cooled toward the glass transition tem- 
perature T g , a drastic slowing occurs for dynamical prop- 
erties such as the structural relaxation time, the diffusion 
constant, and the viscosity, while only small changes are 
detected in static properties [3, HI- Despite the extremely 
widespread use of glass in industry, the formation process 
and dynamic properties of this material are still poorly 
understood. The goal of theoretical investigations of the 
glass transition is to understand the universal mechanism 
that gives rise to the drastic slowing of dynamical prop- 
erties. Numerous studies have attempted to explain the 
fundamental mechanisms of the slowing of the dynamics 
observed in fragile glass (i.e., the sharp increase in vis- 
cosity near the glass transition). However, the physical 
mechanisms behind this slowing have not been success- 
fully identified. 

Recently, "dynamical heterogeneities" in glass-forming 
liquids have attracted much attention. The dynam- 
ics of glass-forming liquids are not only drastically slow 
but also become progressively more heterogeneous upon 
approaching the glass transition. Dynamical hetero- 
geneities have been detected and visualized through sim- 
ulations of soft-sphere systems @-@], hard-sphere sys- 
tems @, and Lennard- Jones (LJ) systems jlQ ]. and 
through experiments performed on colloidal dispersions 
using particle-tracking techniques pdl4l3j . Insight into 
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the mechanisms of dynamical heterogeneities will lead to 
a better understanding of the slowing of the dynamics 
near the glass transition. 

The properties of dynamical heterogeneity can be char- 
acterized by the following three quantities: the correla- 
tion length, the intensity, and the lifetime. In a system 
displaying dynamical heterogeneity, the particles can be 
divided into "slow" and "fast" sub-sets. The slow and 
fast particles form cooperative correlated regions, and 
these slow and fast regions migrate in space over time. 
These three quantities (the correlation length, the inten- 
sity, and the lifetime) can describe the static and dynamic 
properties of the slow and fast regions. The intensity 
measures the average variance of the slow and fast re- 
gions, and the correlation length characterizes the spatial 
extent of the slow and fast regions. The lifetime repre- 
sents the time scale at which the slow and fast regions 
migrate in the space. 

The correlation length, intensity, and lifetime can be 
investigated using the correlation functions of the par- 
ticle dynamics. In fact, we can evaluate the correla- 
tion length £4 and the intensity Xi by calculating the 
four-point correlation functions that correspond to the 
static structure factors of the particle dynamics. Sev- 
eral simulations 

d H ElMiJ experiments HHH, and 
mode-coupling theories |29M31j| have estimated £4 and 
Xi using four-point correlation functions and have re- 
vealed that £4 and Xi increase with decreasing temper- 
ature (or an increase in the volume fraction in the case 
of hard-sphere systems). Furthermore, we can quantify 
the lifetime Ti 10 t or o using the multiple-time extensions of 
the four-point correlation functions (i.e., the multi-time 
correlation functions) that correspond to the time corre- 
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lation functions of the particle dynamics. Recent simula- 
tions have Quantified Thetero using multi-time correlation 
functions [a. l32 - l35| . Various experiments, including pho- 
tobleaching techniques and nuclear magnetic resonance, 
have also measured Thetero [HI l36l - l4l| . It was reported 
that Thetero increases dramatically with decreasing tem- 
perature or an increase in the volume fraction and can 
exceed the a-relaxation time near the glass transition. 

As mentioned above, there have been many studies on 
the correlation length, the intensity, and the lifetime of 
dynamical heterogeneity near the glass transition. How- 
ever, knowledge of and measurements relating to the life- 
time are still limited. Moreover, individual studies have 
been restricted to only the correlation length and the in- 
tensity or only the lifetime, and the relationship between 
the length and time scales of dynamical heterogeneity 
remains controversial despite its importance [HI, [35| . 

The aim of the present study is to examine all the 
three quantities (the correlation length, the intensity, and 
the lifetime of dynamical heterogeneity) consistently. We 
performed molecular-dynamics (MD) simulations and in- 
vestigated dynamical heterogeneity using the correlation 
functions of the particle dynamics. In our previous study 
[42j , we evaluated the lifetime Thetero (t) in two different 
time intervals: the a-relaxation time t q and the time T ngp 
at which the non-Gaussian parameter of the Van Hove 
self-correlation function is maximized. In the present 
study, in addition to the lifetime Thetero (t), we quanti- 
fied the correlation length £i(t) and the intensity X^if) 
from the same order parameter used when calculating 
Thetero (t) . Furthermore, we examined the time-interval 
dependence of the lifetime Thetero (i) to understand the 
lifetime and the dynamic properties of dynamical hetero- 
geneity in more detail. 

The paper is organized as follows. In Sec. |Hl we ex- 
plain the correlation functions of the particle dynamics. 
We show that dynamical heterogeneity can be system- 
atically examined using the correlation functions of the 
particle dynamics. In Sec. IIII1 we briefly review our 
MD simulation and present some results from conven- 
tional density correlation functions. In Sees. IIVI and [V] 
the results for dynamical heterogeneity are presented. In 
Sec. IIVI we first show three quantities, the correlation 
length, the intensity, and the lifetime, which are consis- 
tently calculated from a single order parameter. In Sec. 
IVl we next present the time-interval dependence of the 
lifetime Thetero (t). In Sec. IVI[ we summarize our results. 



II. CORRELATION FUNCTIONS OF PARTICLE 
DYNAMICS 

As we mentioned, dynamical heterogeneity can be 
characterized by three quantities: the correlation length, 
the intensity, and the lifetime. In this section, we intro- 
duce the correlation functions of the particle dynamics 
and demonstrate that these three quantities can be sys- 
tematically evaluated in terms of the correlation func- 
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FIG. 1. (Color) Schematic illustration of the time configu- 
ration of the correlation functions of the particle dynamics: 
(a) the spatial correlation function of the particle dynamics 
in the time interval [0, t] and (b) the time correlation function 
of the particle dynamics between the two time intervals [0, t] 
and [t s +t,t 3 +2t]. 

tions of the particle dynamics. 

The conventional two-point correlation function 
F(k,i) represents the correlation of local fluctuations 
8n(k,t) with some order parameter, such as particle 
density. The expression Sn(k,t) is the Fourier com- 
ponent k of the fluctuations at time t, and F(k,t) — 
(Sn(k,t)Sn(~k,0)), where k = \k\. When t = 0, 
S(k) = F(k,t = 0) is the spatial correlation of Sn(k,0) 
(i.e., the static structure factor), and we can examine the 
static structure of the order parameter by the wavenum- 
ber dependence of S(k). When t > 0, F(k,t) describes 
the particle dynamics in the time interval [0, t], averaged 
over the initial time and space. As the time interval t in- 
creases, F(k, t) decays in the stretched exponential form, 

where r(k) is the relaxation time of the two-point cor- 
relation function that represents the characteristic time 
scale of the average particle dynamics. 

To examine the structure and motion of spatially het- 
erogeneous dynamics, we must calculate the correlation 
of the local fluctuations, SQk(q,to,t), in the particle 
dynamics. The expression SQk(q, to, t) is the Fourier 
component q of the fluctuations in the particle dynam- 
ics that are associated with a microscopic wavenumber 
k in the time interval [to,t + t\. F(k,t) is equal to 
Qk(q,to,t) averaged over the initial time to and space, 
i.e., F(k,t) ~ (Qk(q,to,t)}. The correlation function de- 
fined by 

S 4 ,fc(<7, t) = (SQ k (q, 0, t)5Q k (-q, 0, t)), (2) 

represents the spatial correlation of the particle dy- 
namics in the time interval [0,t]. The time configu- 
ration of Si.k{q,t) is schematically illustrated in Fig. 
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QJa). Si t k(q,t) is the four-point correlation function. 
We can examine the structure of the particle dynamics 
by the wavenumber dependence of S^idq, t). At small 
wavenumbers of q, S^k^q, t) can be approximated by the 
simple Ornstein-Zernike (OZ) form [5l Il6j|. 



III. SIMULATION MODEL AND RESULTS 
FROM THE DENSITY CORRELATION 
FUNCTIONS 



A. Simulation model 



S4 ik (q,t) = 



X4,fe(*) 



q 2 U,k 



■Xt) 1 



(3) 



where £4 fc(t) is the correlation length of the heteroge- 
neous dynamics and Xi,k{t) is the intensity of the het- 
erogeneous dynamics, which is the long- wavelength limit 
of S 4 , k (q,t), i.e., X4,k( t ) = lim <j-K) S i<k (q,t). Note that 
Xi,k(t) is termed the four-poin t dy namical susceptibility. 
According to previous studies |2ll424j . to quantify ^4 lk (t) 
and Xi,k(t) accurately, we need to use a large system 
or determine X&,k(f) through other means (not OZ fit- 
ting). In their work [Ii,[24]], E. Flenner et al. claimed 
that accurate determinations of ^4 lk (t) and Xi,k(t) can 
be made by fitting S4 t k(q, t) to the OZ form in the range 
of q^4,k(t) < 1-5. They used a large system composed of 
8 x 10 4 particles to obtain accurate fits for the OZ form. 
In the present study, we used a large system with 10 5 
particles to fit S4 tk (q,t) to the OZ form accurately and 
obtain the values of £4,fc(£) and X4,k(t)- 

Furthermore, the time correlation function defined by 



F 4 , k {q, t s ,t) = (SQ k (q, t s + t, t)6Q k (-q, 0, t)>, (4) 



represents the correlation of the particle dynamics be- 
tween the two time intervals [0, t] and [t a + t,t s + 2t]. 
The value t s is the time separation between the two time 
intervals [0,i] and [t s + t,t 3 + 2t]. The time configu- 
ration of i^4,fc(g, t s , t) is schematically illustrated in Fig. 
[lib). F4,k(q,t s ,t) is the multi ple- time extension of the 
four-point correlation function [34l . l35j . As the time sep- 
aration t a increases, F 4 ^ k (q,t s ,t) with fixed t decays in 
the stretched exponential form, 



F^ k {q,t s ,t) 
F 4 , k (q,0,t) 



exp 



TA,k{q,t) 



(5) 



where T 4 k (q,t) is the relaxation time of the correlation 
of the particle dynamics. We determined the lifetime 
Thotoro(i) of the heterogeneous dynamics as T4. k (q,t) at 
q = 0.38. We need longer trajectories of the simulations 
to quantify Th c tero(£) than to quantify £,4, k (t) and X4,fc(t)- 
In the present study, we used a smaller system with 10 4 
particles to calculate Thetero(i)- As explained above, we 
can systematically evaluate the correlation length, the 
intensity, and the lifetime of dynamical heterogeneity by 
calculating the correlation functions of the particle dy- 



We performed MD simulations in three dimensions on 
binary mixtures of two different atomic species, 1 and 
2, with a cube of constant volume V as the basic cell, 
surrounded by periodic boundary image cells. The par- 
ticles interacted via their soft-sphere potentials, v a b(r) = 
e(cr a b/r) 12 ; where r is the distance between two particles, 
Cab = (c a + fb)/2, and a, b £ 1,2. The interaction was 
truncated at r = 3a a b- The mass ratio was m%lm\ = 2, 
and the diameter ratio was 02/01 = 1-2. This diameter 
ratio avoided system crystallization and ensured that an 
amorphous supercooled state formed at low temperatures 
[43l | . As mentioned in Sec. [Til we used two systems: a 
small system with Ni = N 2 = 5 x 10 3 (N = Ni + N 2 = 
10 4 ) particles and a large system with Ni — N 2 — 5 x 10 4 
(JV = ^1 + ^2 = 10 5 ) particles. The large system was 
used to quantify the correlation length £4,fe(t) and the in- 
tensity Xi,k(t), and the small system was used to quantify 
the lifetime Thotero(t)- In the present paper, the following 
dimensionless units are used: length, a\\ temperature, 
e/fcg; and time, t = (m\o\i 'e) 1 / 2 . The particle density 
was fixed at the high value of p = N/V = 0.8. The sys- 
tem lengths were L = V 1 / 3 = 23.2 and 50.0 for the small 
and large systems, respectively. Simulations were per- 
formed at T = 0.772, 0.473, 0.352, 0.306, 0.289, 0.267, 
and 0.253. Note that the freezing point of the corre- 
sponding one-component model is approximately T = 
0.772 (r off = 1.15) HI]. Here, r eff is the effective density, 
which is a single parameter characterizing this model. At 
T = 0.253 (r c ff = 1.52), the system is in a highly super- 
cooled state. We used the leapfrog algorithm with time 
steps of 0.005 when integrating the Newtonian equation 
of motion. At each temperature, the system was care- 
fully equilibrated under the canonical condition so that 
no appreciable aging effect was detected for various quan- 
tities, including the pressure and the density correlation 
function. Once equilibrium was established, data were 
taken under the microcanonical condition. The length 
of the data collection runs was at least 100 times the a- 
relaxation time, T a , for the small system and 10 times r a 
for the large system. Information regarding this model, 
such as the static structure factor, the intermediate scat- 
tering function, and the mean square displacement, can 
be found in previous works [H, [44| . 



B. Single- particle and collective-particle diffusive 
motion 



Before showing the results for dynamical heterogene- 
ity in supercooled liquids, we present some results from 
an investigation of the average particle dynamics us- 
ing conventional density correlation functions. Let us 
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consider the density correlation functions F sa (k,t) and 
F a (k,t) (a e 1,2), denned by 



(6) 



F„(M) = <<MM)<M-M)>, 



where 5n a j{k,t) = exp[— ik • r a j(t)] is the Fourier com- 
ponent k of the tagged particle density fluctuations of 
particle species a, and dn a (k, t) = ^2f=i exp[— zfc • r a j(t)] 
is the Fourier component k of the density fluctuations of 
particle species a. The terms F sa (k,t) and F a (k,t) de- 
scribe the single-particle and collective-particle motion, 
respectively [45|]. We calculated F sa (k,t) and F a (k,t) for 
a wide range of wavenumbers, k — 0.35 — 40. As seen 
in Eq. ((T|), F sa (k,t) and F a (k,t) decay in the stretched 
exponential form, 



F 3 g(k,t) 

F sa (k,0) 

F a (k,t) 
F a (k,0) ' 



exp 



exp 




(7) 



where T sa (k) and r ca (k) are the wavenumber-dependent 
relaxation times of F sa (k,t) and F a (k,t), respectively. 

In Fig. [21 T sa {k) and r ca (k) are plotted for parti- 
cle species 1 and 2 as functions of the wavenumber k. 
In Fig. [UJa), T sa (k) approaches r sa (k) = D~ c }k~ 2 at 
small wavenumbers fc, where D sa is the diffusion constant 
of the single-particle motion of particle species a. The 
term D sa is calculated by D sa = lim t ^ 00 ([Ar a (i)] 2 )/6i; 
where ([Ar a (t)] 2 ) is the mean square displacement of 
particle species a; ([Ar a (t)} 2 ) = (Zf=il Ar aj{t)} 2 /N a ); 
Ar a j(t) — r a j(t) — r a j(0). The diffusion constant of par- 
ticle species 1 is larger than that of particle species 2 
at every temperature. However, T ca {k) also approaches 
T ca (k) ~ fc~ 2 at small k in Fig. [2fb). This behavior indi- 
cates that the collective-particle motion is also diffusive 
at large length scales, and we can define the diffusion 
constant of the collective-particle motion D a as 



D„ 



lim T rn 1 (k)k 



We can see that F)\ and Z?2 are almost identical between 
D sl and D s2 in Fig. 02b). 



C. Correlation between the static structure factor 
and the relaxation time 

Next, we consider a density variable defined by 

p eB (k, t) = ofrai(k, t) + cr|n 2 (fe, t), (9) 

that represents the degree of particle packing and is the 
effective one-component density of our binary mixture 
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FIG. 2. (Color) The wavenumber dependence of (a) T 3a (k) 
and (b) r ca (k) for particle species a = 1 and 2. Temperatures 
are 0.772 — 0.253 from the lowest curve to the highest. The 
black line is T sa (k) — D^fe -2 , where D sa is the diffusion con- 
stant of the single-particle motion calculated from the mean 
square displacement. 



[5j. We calculated the static structure factor S pp (k) and 
the time correlation function F pp (k,t), defined by 



Sp P {k) = ^{5p cS (k,0)5p cS (-k,0)) : 
F pp {k,t) = (S PcS {k,t)Sp ctf {~k,0)). 



(10) 



The terms S pp (k) and F pp {k, t) represent the static struc- 
ture and the average particle dynamics of the effective 
(8) one-component fluids, respectively. Note that F pp (k,t) 
decays in the stretched exponential form with the relax- 
ation time T pp (k), as in Eq. ([T]). 

Figure [3] shows the wavenumber dependence of S pp (k) 
in (a) and T pp (k) in (b). We can see that S pp (k) and 
Tpp(fc) are maximized and minimized at almost the same 
wavenumbers k. For example, both S pp (k) and T pp (k) 
have first-peak values around k = lix. Thus, there is a 
correlation between S pp (k) and r pp (fc), i.e., between the 
static structure and the particle dynamics. A similar 
correlation was found in the Lennard- Jones (LJ) model 
(irjj . wat er |47| . and the polymer model According 
to Ref. 




(fc) is modulated by S pp (k), and this 



modulation can be understood as a consequence of 
Gennes narrowing" [49(. 
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FIG. 3. The wavenumber dependence of (a) S pp (k) and (b) 
Tp P (k). The temperatures are 0.772 — 0.253 from the lowest 
curve to the highest. 




FIG. 4. Time intervals r a and T ngp versus the inverse tem- 
perature 1/T. We use these two time intervals to define the 
local dynamics. 



IV. RESULTS I: CORRELATION LENGTH, 
INTENSITY, AND LIFETIME 

In this section, we present the results of three quan- 
tities: the correlation length, the intensity, and the life- 
time of dynamical heterogeneity. We consistently evalu- 
ated these three quantities using a single order parame- 
ter representing the particle dynamics and its correlation 



functions. In our previous study [42[, we evaluated the 
lifetime ih otcro (i). For this study, in addition to the life- 
time Thotcro(i); we quantified the correlation length £i(t) 
and the intensity X4(i) from the same order parameter 
used when calculating Thctcro(i)- As in our previous study 
we used two time intervals, r Q and r ngp , to define 
the local dynamics. The value r Q is the a-relaxation 
time defined by F s ±(k m ,T a ) = e -1 , where F s i(k,t) is the 
self-part of the density time correlation function for par- 
ticle species 1 as defined in Eq. ^ and k m — 2n is the 
first-peak wavenumber of the static structure factor. The 
value T ngp is the time at which the non-Gaussian param- 
eter aa(t) of the Van Hove self-correlation function 
defined as a 2 (t) = 3([Ari(t)] 4 )/5([Ari(i)] 2 ) 2 -l is maxi- 
mized. In Fig. 21 we show r Q and T ngp as functions of the 
inverse temperature 1/T. At T = 0.306, t q ~ T ngp , and 
T a grows exponentially larger than r ngp with decreasing 
temperature at T < 0.306. This tren d ag rees with other 
simulation results for L J systems [H], [52| . 



A. The structure of the heterogeneous dynamics 

We first examined the structure of the heterogeneous 
dynamics. We calculated the displacement of each 
particle of species 1 in the time interval [to, to + t]; 
Arij(to,t) =ry(«d+t)-r y (to) (j = l,2,...,iVi), and 
the particle mobility (to , t) of each particle was defined 
as 



[Ary(to,t)] s 



(11) 



<[Ary(t 0! i)] 2 )- 

In Fig. [SJ we show the spatial distribution of the particle 
mobility at T — 0.253, which is the lowest temperature 
in our simulations. In the figure, the particles are drawn 
as spheres with radii a 2 ^ (t Q , t) located at 



Rij{to,t) = \ [rij{to) 



■n&o + t)]. 



(12) 



Notice that a\j{to,t) > 1 (a\j(to,t) < 1) means that 
the particle j moves more (less) than the mean value of 
the single-particle displacement, i.e., particle j is mobile 
(immobile). In Fig. [5j the red (blue) spheres represent 
ajitoit) > 1 (a|(io,£) < 1). In the figure, the two hetero- 
geneity structures in r Q and r ngp are both significant but 
differ considerably. There are many more red spheres in 
[SJa) than in[SJb), but there are much larger red spheres 
in [5j^b) than in[5ja). This pattern means that many mo- 
bile particles contribute to the heterogeneity in the time 
interval T a , but in the time interval r ngp , relatively few 
particles are mobile and contribute to the heterogeneity. 

We also calculated the spatial correlation function of 
the particle dynamics expressed in Eq. @. We con- 
sidered the local fluctuations in the particle mobility of 
particle species 1 defined by 



3 = 1 



Rij(t ,t)), (13) 
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(a) t = r a 



(b) * = T ngp 





FIG. 5. (Color) The distribution of the particle mobility 
a\j{to,t) for particle species 1. The time intervals are (a) 
[t , t + r a ](t = r a ) and (b) [to, to + T ngp ](i = r ngp ). The 
temperature is 0.253. The radii of the spheres are alj(to,t), 
and the centers are at Rij(to,t). The red and blue spheres 
represent a\j(to,t) > 1 (mobile particles) and a\j{to,t) < 1 
(immobile particles), respectively. 
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ter 5Vi(q,to,t) represents the local fluctuations in the 
particle dynamics in the time interval [to, to + t]. Note 
that the same order parameter as that of 5T>i(q, to, t) is 
used to evaluate the lifetime in our previous study |42|. 
We used ST>i(q,t ,t) as SQk{q,t ,t) in Eq. ([2]), and the 
correlation function defined by 

Si(q,t) = ^r(S7h(q,Q,t)S1h(-qAt)), (15) 
iVi 

corresponds to Si,k{c[,t). The term S4,(q,t) represents 
the spatial correlation of the particle dynamics in the 
time interval [0,t]. The time configuration of S±{q,t) 
is schematically illustrated in Fig. [lja). We were able 
to examine the structure of the heterogeneous dynamics 
using the wavenumber dependence of S±{q,t). 

Figure |5] shows the wavenumber dependence of S^{q, t) 
of particle species 1 for t = r a and T ngp . We cal- 
culated Si(q,t) using a larger system with N = 10 5 
to quantify the correlation length and the intensity ac- 
curately. At small wavenumbers of q (long-distance 
scales) , the correlations in r Q and r ngp both become large 
with decreasing temperature in a similar manner. How- 
ever, at large wavenumbers of q (short-distance scales), 
54(g, r ng p) grows larger than Si(q,T a ) at low tempera- 
tures, which reflects that more highly mobile particles 
exist in the time interval be seen in the vi- 

sualization shown in Fig. [5] 
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FIG. 6. The wavenumber dependence of S4(q,t) for particle 
species 1 at T — 0.772 — 0.253. The time interval t is r Q in 
(a) and T ngP in (b). Note that S4,(q,t) was calculated with a 
larger system, TV = 10 5 . 



or the Fourier component q of 6T>i(r,to,t), 

ATi 

5V 1 (q,to,t) = J2 Sa lj( t o,t)eM-iq-Rij{h,t)], (14) 
3=1 

where &a\j(t ,t) = a\j(to,t) — 1 is the fluctuation of the 
particle mobility of the particle j. The order parame- 



B. The correlation length, the intensity, and the 
lifetime of the heterogeneous dynamics 

We next quantified the correlation length and the in- 
tensity of the heterogeneous dynamics in t q and T ngp . 
As mentioned for Eq. ([3]), at small wavenumbers of q, 
S4,(q,t) can be approximated by the simple OZ form, 



S 4 (q,t) 



(16) 



where £i(t) is the correlation length and Xiif) is the in- 
tensity. The values £n(t) and Xi(f) correspond to £,i,k(t) 
and X4,fc(i) in Eq. (|3]), respectively. To obtain accurate 
values of £i(i) and Xi{t)i Si{q,t) was carefully fitted to 
the OZ form in the range of q£,±(t) < 1.5 [23I. l24j|. 

Figure [7] shows £4(4) versus t q in[7Ja) and Xi{t) versus 
£,i(t) in[7Jb) for t = r a and T ngp . For the time interval 
t = r a , we examined the scaling relationships between 
T a , £i(T a ), and Xi{ T a)- As in Fig. [TJa), we obtained a 
power law fit £4(t q ) ~ tU z with l/z = 0.1. The scaling 
exponent 0.1 is very small. We also found that a fit to 
£,4,(tc) ~ (liiTa) 1 / 61 with 1/9 = 0.77 provides a better de- 
scription of the data over a larger range of T a , so there 
is a slower-than-power law increase of £4(T a ) with r a . 
Therefore, ^{T a ) increases much more slowly compared 
to Tq as the temperature decreases. This result is consis- 
tent with the most recent results 0, [24[ ■ Furthermore, 
we obtained the scaling relationship X4(t q ) ~ ^4(T a ) 2 ~ v 
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FIG. 7. (a) The correlation length versus T a . (b) The 
intensity x-i{p) versus the correlation length (,4,(t). The time 
intervals are t = t q and T ngp . The straight lines are power 
law fits: £ 4 (t q ) ~ r°' 1±0 - 01 and X4 (t„) ~ £4(T Q ) 3 ' 2±ai . The 
dashed curve is a fit to ^4(r Q ) ~ (lnr a )°' 77±0 12 . 



with 2 — i] = 3.2 in Fig. [7|b). Our scaling exponent 3.2 
agrees well with Ref . [23| . 

On the other hand, when the time interval t is T ngp , 
^4(T„gp) reaches a value near r a ~ 250 and then plateaus 
as t q increases in Fig. Eta), while X4(r ng p) monotoni- 
cally gets large in Fig. [T^b). Here, we remark that even 
in the range where ^4(r ng p) plateaus, X4(r ngp ) contin- 
ues to increase with r Q . This behavior of ^4(T ngp ) and 
X4( r n g p) m ay be due to the choice of the order parame- 
ter ST>i(r,to,t). ST>i(r,to,t) represents the distribution 
of the square particle displacements, so the mobile parti- 
cles with large displacements contribute to the structure 
factor S±(q,t) much more than the immobile particles. 
At t = r ngp and at low temperatures, our simulation 
data show that only a small number of particles tend 
to have very large displacements compared to the other 
mostly immobile particles, which can be also seen in Fig. 
O In this situation, the intensity X4( T n gP ) can increase 
from the contribution of particles with very large dis- 
placements, whereas the correlation length ^4(r ngp ) can- 
not increase because the number of mobile particles in- 
volved in the correlated regions decreases. 

So, we used an another order parameter SV\ og i (r, to, t) 



FIG. 8. The correlation length £t(t) and the intensity x^if) 
calculated by using SV\ ogl (r,to,t), i.e., the distribution of 
the logarithm of the square particle displacements, (a) £4(t) 
versus T a . (b) Xi(t) versus £i(t). The time intervals are t — r a 
and Tngp. The straight lines are power law fits: ^4(r Q ) ~ 
t°- 1±0 ' 01 and X 4(Vc0 ~ £4(r Q ) 3 - 5±0 ' 2 . The dashed curve is a 
fit toe4(r a )~(lnr Q )°- 91±0 - 12 . 



defined as 

8V log i (r, t , t) = ^ln a^to, t)S(r - J*y(*o, *))■ (17) 

The order parameter 5T>i og i(r, to, i) represents the distri- 
bution of the logarithm of the square particle displace- 
ments. We quantified and Xi(t) m terms of the 
structure factor of 5T>\ og i(r,to,t). Notice that in this 
case, the immobile particles with small displacements 
contribute to the structure factor to the same extent 
as the mobile particles because of the logarithm oper- 
ation. Figure [5] shows ^(t) and Xi{t) calculated by 
i og i(r, to, i). At the time interval t = r a , the scaling 
relationships between r Q , ^4(r Q ) and X4( T a) are almost 
same as those examined by 5V\(r,to,t). On the other 
hand, at the time interval t = T ngp , both ^4(r ngp ) and 
X4( r n gP ) reach values near r Q ~ 250 and then plateau 
as r a increases, i.e., both ^4(r ngp ) and X4( r n gP ) show 
the plateau of the heterogeneity. This result indicates 
that the behavior that X4( r n gP ) continues to increase 
in Fig. [TJa) is due to the choice of the order parame- 
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ter ST>x(r,to,t). For the confirmation purpose, we also 
calculated the lifetime Thctcro (i) by using ST>i og i(r,to,t) 
and checked that the lifetimes of both ST>i(r,to,t) and 
ST> i gi(f*, to, t) have almost same values and behave in 
the same manner with decreasing temperature. 

In our previous study, we already determined the scal- 
ing relationships between the lifetime Thctcro (i) and r a , 
i.e., Thetero(Ta) ~ t^ 08 and T hetoro (r ngp ) ~ r°- 91 . At the 
time interval t = r a , we can determine the scaling rela- 
tionship between the correlation length and the lifetime: 

Thctcro (T a )~Ur a ) W - 8 , (18) 

or 

Thctcro(T Q ) - exp(fc^ 4 (r Q ) 1 - 3 ). (19) 

The scaling exponent 10.8 of the power law scaling is very 
large, and there is an exponential growth of Thctcro (t q ) 
with ^4(r Q ). Furthermore, at the time interval t — T ngp , 
although the correlation length £4(r ngp ) plateaus, the life- 
time Thetero(T"ngp) continues to increase dramatically with 
decreasing temperature. Thus, we can conclude that 
the lifetime Thctcro gets large dramatically with decreas- 
ing temperature, whereas the correlation length £4 and 
the intensity Xi increase slowly compared to Thctcro or 
plateau, i.e., the time scale of dynamical heterogeneity 
grows faster than the length scale and the intensity of 
dynamical heterogeneity. 



(a) t = O.OOItc (b) t = 0.1r ( 




FIG. 9. (Color) The time-interval dependence of the distri- 
bution of the particle mobility a\j{to,t) for particle species 
1. The time intervals are (a) t — 0.001r Q , (b) t — 0.1r a , 
(c) t — r a and (d) t = 100r a . The temperature is 0.267. 
The radii of the spheres are a\j(ta,t), and the centers are at 
Rij(to,t). See also the caption for Fig. [5] 



V. RESULTS II: TIME-INTERVAL 
DEPENDENCE OF LIFETIME 

In this section, we present the results for the time- 
interval dependence of the lifetime Thctcro (t)- To exam- 
ine the lifetime in more detail, we evaluated the lifetime 

Thctcro 

(t) for various time intervals t and determine how 
the lifetime Thctcro (t) depends on the time interval. 

A. The structure of the heterogeneous dynamics 

First, we examined the time-interval dependence of the 
structure of the heterogeneous dynamics. The spatial 
distribution of the particle mobility a^(io,i) for particle 
species 1 is shown in Fig. [5J in which the time interval 
increases from t = 0.001t q to IOOtq,. The temperature is 
0.267. We can recognize that the heterogeneity is much 
weaker at the short time interval t = 0.001t q in[9]Ja). As 
the time interval t increases, the heterogeneity increases 
and is maximized between t = 0.lT a in |9jb) and t — 
t q in[9]Jc). As the time interval t increases further, the 
heterogeneity decreases and is weakened at t = 100T a in 
Hd). 

We also calculated the spatial correlation function 
Si(q,t) defined in Eq. (fT5]l . Figure [10] shows the 
wavenumber dependence of S^q, t) for various time inter- 
vals t. We can see that as the time interval t gets large, 



54 (g, t) increases, is maximized, and then decreases in 
the region of q. This behavior agrees with the visualiza- 
tion of the heterogeneity structure shown in Fig. [9j It is 
also seen that at large time intervals, S^q, t) maintains a 
constant value, independent of the wavenumber g. This 
result indicates that at large time intervals, particle mo- 
bilities are uniformly distributed throughout space, i.e., 
the structure of the particle dynamics is spatially homo- 
geneous. Note that the static structure factor 5n(g) also 
becomes constant at small g (long-distance scale). 



B. The motion of the heterogeneous dynamics 

Next, the motion of heterogeneous dynamics was in- 
vestigated in detail. In our previous study (42|, we in- 
vestigated the motion of the heterogeneous dynamics at 
the time intervals t = r a and T ngp . Our results suggested 
that the heterogeneous dynamics might migrate in space 
with a diffusion-like mechanism. In the present study, we 
examined the motion of the heterogeneous dynamics at 
various time intervals t. As in our previous study [42j . 
we calculated the time correlation function of the parti- 
cle dynamics expressed in Eq. (@|. We used 5T>i(q,t Q ,t) 
as SQk(q,to,t) in Eq. ((4]), and the correlation function 
S v {q,t s ,t), 

Sv(q,t s ,t) = («Ji(g,t a + t,t)«M-g,0,t)> ) (20) 
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FIG. 10. The wavenumber dependence of S&{q,t) at T = 
0.306. The time intervals are (a) 0.05r a , 0.lT a , 0.5r a , and r a 
from the lowest curve to the highest, and (b) T a , 5r a , 10T a , 
50r a , 100r Q , and 300r a from the highest curve to the lowest. 
Si(q,t) is maximized at t = r a . The dashed curve represents 
the static structure factor Sn(q). 



corresponds to F4 t k(q,t s ,t) and represents the correla- 
tion of the particle dynamics between two time inter- 
vals [0,t] and [t s + t,t s + 2i\. The time configuration of 
Sx>(q,t s ,t) is schematically illustrated in Fig. [ljb). As 
mentioned for Eq. ([5]), when the time separation t s in- 
creases, S-p(q,t s ,t) with fixed t decays in the stretched 
exponential form, 



Sy(q,t Sl t) 
Sv(q,0,t) 



exp 



Th(q,t) 



(21) 



where r^(g, t) is the wavenumber-dependent relaxation 
time of Sj)(q, t s , t) and represents the time scale at which 
the heterogeneous dynamics move in the space. The ex- 
pression Th{q,t) corresponds to T^h{q,t) in Eq. §5§. We 
were able to examine the motion of the heterogeneous 
dynamics using the wavenumber dependence of Th (q, t) . 

In Fig. Ill) we show the wavenumber dependence of 
Th{q,t) for various time intervals t. As the time inter- 
val t increases, at small wavenumbers q, Th{q, t) increases 
monotonically and approaches the relaxation time T c i(q) 
of the two-point density correlation function shown in 
Fig. [2] This result indicates that the heterogeneous dy- 
namics at large time intervals behave like the particle 



FIG. 11. The wavenumber dependence of i~h(q,t) at T = 
0.306. The time intervals are 0.05r a , 0.1r Q , 0.5r a , T a , 5r a , 
10r a , 30t q , 50r a , 70r a , 100r a , and 300r Q from the lowest 
curve to the highest. The dashed curve is the relaxation time 
Tci (q) of the two-point density correlation function. 




FIG. 12. The wavenumber dependence of Th s {q,t) at T = 
0.306. The time intervals are 0.05r a , 0.1r Q , 0.5r a , r a , 5r a , 
IOtc, 30t q , 50r a , 70r Q , 100r Q , and 300r Q from the lowest 
curve to the highest. The dashed curve is the relaxation 
time Tsi(q) of the self-part of the two-point density correlation 
function. 



density, the motion of which is diffusive at a long dis- 
tance scale. Notice that Th{q : t) at large t and small q is 
proportional to g -2 , which indicates this diffusion mech- 
anism. We also checked that Sx>(q,t s ,t) decays in the 
exponential form at large t and small q, as does the two- 
point density correlation function. 

Furthermore, we calculated the self-part of Sx>(q, t s ,t), 
defined by: 

Sua (?,*.,*) - (^^ l SVij{q,t.+t,t)SV lj (q,0,t) ), 



where 



(22) 



SVi j (q,h,t)=5c^to,t)exp[-iq-Ri j (t ,t)]. (23) 

The correlation function Sj) s (q,t s ,t) represents the cor- 
relation of the individual particle dynamics between two 
time intervals [0, t] and [t s +t,t s + 2i\. Like Su(q,t s ,t), 
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as the time separation t s increases, Sr> s (q, t s , t) decays 
in the stretched exponential form with the relaxation 
time Th s (q,i). The value of r^ s (g,i) is the time scale 
at which the individual heterogeneous dynamics move in 
space. Figure [12] shows the wavenumber dependence of 
T~hs(Qit) for various time intervals. As the time inter- 
val t increases, Th s (q,t) increases monotonically at small 
q and approaches the relaxation time T s i(q) of the self- 
part of the two-point density correlation function shown 
in Fig. [2] Like Th(q,t), Th s (q,t) is also proportional to 
q~ 2 at large t and small q. Thus, both the collective- 
particle behavior and the single-particle behavior of the 
heterogeneous dynamics at large time intervals are diffu- 
sive, like the motion of particle density at long-distance 
scales. 

That Th(q,t) and T} ls (q,t) approach the relaxation 
times T c \{q) and T s \(q) of the two-point density corre- 
lation functions can possibly be interpreted as follows. 
The correlation function Sx> (<7, t s , t) can be written as 

S v (q, t a , t) = / Sa h (*«+*> l ) 6a lk (0, t) x 

X j=i fc=i (24) 

exp[-iq • (R xj {t s +t,t) - R lk (0,t))} 

From this expression, there can be two types of relaxation 
of Sx>(q,t s ,t) with increasing time of separation t s . One 
is the relaxation due to fluctuations in the particle mo- 
bility, i.e., due to the term «5a\ j {t a +t,t)5a\ k {Q,t) n . The 
other is the relaxation due to particle motion, i.e., due to 
the term "exp[— iq- (Rij(t s +t,t) — i?ifc(0, £))]" , which is 
the same as the relaxation of the density correlation func- 
tions. So, the total relaxation time Th(q,t) of Sj}(q,t s ,t) 
is determined by the two times Ts a (t) and T&(q), where 
Ts a (t) is the time scale at which particle mobility fluctu- 
ates, and r c i(q) is the time scale of particle motion and 
is the relaxation time of the density correlation function. 
Here, for simplicity, we assume that S-z>{q,t s ,t) relaxes 
with these two times T$ a (t) and T c i(q) in the form 



Sy{q,t s ,t) 
Sv(q,0,t) 



exp 



TSa(t) 



x exp 



Td(q) 



The relaxation time Th(q,t) is then mainly determined 
by the smaller time scale of rs a (t) and T c \{q). When 
the time interval t is small, particle mobilities fluctuate 
faster than the particles move, which means that Tg a (t) is 
smaller than r cl (q). In this case, T^(g, t) is determined by 
the time Tg a (t). As the time interval t increases, then the 
fluctuations in particle mobilities becomes slower, and 
the lifetime {q, t) increases accordingly. As the time 
interval t increases further, particle mobilities fluctuate 
more slowly than the particles move, which means that 
T5a{t) becomes larger than r cl (q). In this case, Th(q,t) 
is determined by the time r c i(g) instead of Tg a (t)- Thus, 
at large time intervals, because particle mobilities fluctu- 
ate very slowly, the relaxation time Th(q, t) can be deter- 
mined by the time scale of particle motions, which is the 




KT 3 10 10 10 10 10 10' 
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FIG. 13. The time-interval dependence of Thotero(i). Temper- 
atures are 0.473 in (a), 0.306 in (b), 0.267 in (c), and 0.253 
in (d). The time interval t and the lifetime Th e tero(£) are nor- 
malized by T a - The dotted line indicates the value of r c i(q) at 
q = 0.38, where r c i(q) is the relaxation time of the two-point 
density correlation function. 



relaxation time of the density correlation function. 



C. The lifetime of the heterogeneous dynamics 

Finally, we determined the lifetime Thetero(i) of the het- 
erogeneous dynamics as Th(q,t) at q = 0.38. We show 
the time-interval dependence of the lifetime Th e tero(i) f° r 
various temperatures in Fig. [13] As the time interval t 
increases, T^ c t cm (t) increases monotonically, and at large 
time intervals, Thetero(i) approaches and is limited to the 
relaxation time T c \(q = 0.38) of the two-point density 
correlation function. At large time intervals for which 
Thctcro(£) plateaus, the heterogeneous dynamics migrate 
in space with a diffusion mechanism like that of particle 
density, as we showed in Figs, (jlip and (fT2| . 



VI. SUMMARY 

In this study, we have investigated three quantities 
that characterize dynamical heterogeneity: the corre- 
lation length £i(t), the intensity X4(£), and the life- 
time Thetero(i)- The intensity Xi(t) measures the aver- 
age variance of the slow and fast regions, and the cor- 
relation length £i(t) characterizes the spatial extent of 
the slow and fast regions. The lifetime Th c tcro(i) rep- 
resents the time scale at which the slow and fast re- 
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gions migrate in space. We evaluated all three quan- 
tities using a single order parameter representing the 
particle dynamics and its correlation functions. To de- 
fine the particle dynamics, we used two time intervals, 
t = r a and T ngp . We found that at low temperatures, the 
lifetime Thetero(i) increases dramatically with decreasing 
temperature. In contrast, the correlation length £<i(t) 
and the intensity Xi(t) increase slowly compared to the 
lifetime or plateau. At the time interval t = r a , we 
obtained the scaling relationships between ^(t), Xi{t), 
and T hotoro (T Q ): Xi{T a ) ~ U( T a) 3 ' 2 and r hctcro (r a ) ~ 
£4(Va) la8 or T hctC ro(T Q ) ~ exp(fc£ 4 (r a ) 13 ) . The scaling 
exponent 10.8 of T h ctcro(T a ) ~ £4(t q ) 10 ' 8 is very large, 
and there is an exponential growth of Thctcro (t q ) with 
^(Ta). Furthermore, at the time interval t = T ngp , al- 
though the correlation length ^4(r ng p) plateaus as the 
temperature decreases, the lifetime Thctcro(T ng p) contin- 
ues to increase dramatically. Thus, we can conclude that 
the lifetime Thctcro gets large dramatically with decreas- 
ing temperature, whereas the correlation length £4 and 
the intensity \a increase slowly compared to Thctcro or 
plateau, i.e., the time scale of dynamical heterogeneity 
grows faster than the length scale and the intensity of 



dynamical heterogeneity. 

Furthermore, we investigated the time-interval de- 
pendence of the lifetime Thctcro (t). As the time inter- 
val t increases, Thctcro (t) increases monotonically. At 
large time intervals, the lifetime Thctcro (t) approaches 
and is limited to the relaxation time of the two-point 
density correlation function. At those large time in- 
tervals, the wavenumber-dependent lifetimes Th(q,t) at 
small wavenumbers q (long-distance scales) almost co- 
incide with the relaxation time T c i(g) of the two-point 
density correlation function and are proportional to q~ 2 . 
Therefore, the heterogeneous dynamics migrate in space 
with a diffusion mechanism like that of particle den- 
sity. Note that at large time intervals, particle mobili- 
ties afj (to , t) are uniformly distributed in space, and the 
heterogeneity structure is much weaker. 
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